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In a recent contribution [arXiv:0904:4151] entanglement renormalization was generalized to 
fermionic lattice systems in two spatial dimensions. Entanglement renormalization is a real-space 
coarse-graining transformation for lattice systems that produces a variational ansatz, the multi-scale 
entanglement renormalization ansatz (MERA), for the ground states of local Hamiltonians. In this 
paper we describe in detail the fermionic version of the MERA formalism and algorithm. Starting 
from the bosonic MERA, which can be regarded both as a quantum circuit or in relation to a coarse- 
graining transformation, we indicate how the scheme needs to be modified to simulate fermions. To 
confirm the validity of the approach, we present benchmark results for free and interacting fermions 
on a square lattice with sizes between 6x6 and 162 x 162 and with periodic boundary conditions. 
The present formulation of the approach applies to generic tensor network algorithms. 
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I. INTRODUCTION 

The simulation of strongly correlated fermions in two 
dimensions remains one of the biggest challenges in com- 
putational physics. Quantum Monte Carlo is very pow- 
erful in solving (unfrustrated) bosonic problems, but it 
fails for fermionic systems because of the negative sign 
problem^ which implies an exponential scaling of the 
computational cost with system size and inverse temper- 
ature. Accurate simulations of fermions are crucial to 
gain further insight into phenomena where strong corre- 
lations play an essential role, such as high-temperature 
superconductivity and the fractional quantum Hall effect. 
Progress in this direction has been made in recent years 
with various other methods such as the cluster dynamical 
mean-field theory^ variational Monte CarloP Gaussian 
Monte CarloJ^ and diagrammatic Monte CarloP How- 
ever, even the phase diagram of the simplest lattice model 
of strongly correlated electrons, the Hubbard modelPis 
still controversial. 

One dimensional fermionic problems can be accurately 
solved by the successful density matrix renormaliza- 
tion group (DMRG) method. 13 However, DMRG-type ap- 
proaches fail for large systems in two dimensions be- 
cause of an accumulation of short-range entanglement 
across block boundaries under successive renormaliza- 
tion group (RG) transformations. In recent years several 
ideas to extend DMRG to higher dim ensions by means of 
tensor networks have been developed P 9 | 10 | 11 | 12 | 13 | 14 | 15 | :L6 1 
We focus here on one particular class of tensor net- 
works called the multi-scale entanglement renormaliza- 
tion ansatz (MERA), which is based on the concept of 
entanglement renormalization^ The key idea is to ap- 
ply unitary transformations (disentanglers) locally to the 
system in order to remove short-range entanglement be- 
fore each coarse-graining step. This prevents the accumu- 
lation of degrees of freedom under successive RG trans- 
formations, so that arbitrarily large lattice sizes for crit- 
ical and non-critical systems in one and two dimensions 
can be addressed. The MERA is a variational ansatz of 



the ground-state (or low energy subspace) of a system, 
from which arbitrary local observables and two-point cor- 
relators can be easily extracted. The accuracy of the 
ansatz depends on the amount of entanglement in the 
system, and can be controlled by a refinement parame- 
ter x- I n one-dimensional lattices, the sche me has b een 
used to study several quantum spin systems ^ 7 * 18 * 19 * and 
shown to be particularly suited to study quantum crit- 
ical pointg^ 17 * 19 * 20 * 2 ^ 2 ^. In two dimensions, accurate 
results have b een obtained for free fermionic a nd boso nic 
systems) 21 * 22 ! as well as quantum spin systems) 23 * 24 * 25 ! in- 
cluding large lattices beyond the reach of exact diagonal- 
ization and DMRGpl and frustrated antiferromagnets 
beyond the reach of quantum Monte CarloP^ In addi- 
tion, an analytical MERA characterization has been pro- 
vided for the gro und st ates of a large class of models with 
topological orderP3221 

In a recent paper J2HI entanglement renormalization and 
the MERA were generalized to fermionic systems. Here 
we present a more detailed description of the fermionic 
MERA and provide additional benchmarking results for 
free and interacting fermions in two dimensional lattices. 
The paper is organized as follows. In Sec. |ll]we overview 
the MERA formalism as a means to prepare its gener- 
alization to fermionic systems. The MERA is presented 
both as a quantum circuit and as implementing a coarse- 
graining transformation. Practical calculations involve 
contracting diagrams or tensor networks. These corre- 
spond to different elements (such as the ascending and 
descending superoperators and environments) needed in 
order to compute expectation values from the MERA or 
to optimize this variational ansatz. 

In Sec. |III| we introduce the two incredients necessary 
in order to represent and simulate fermions. First, the 
tensors that constitute the MERA, namely disentanglers 
and isometries, are chosen to be parity invariant or Z 2 
symmetric. Z 2 symmetric tensors are convenient in order 
to account for parity preservation. Second, we associate 
a fermionic swap gate to every crossing of lines in a di- 
agram. This gate accounts for fermionic statistics, and 
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is the key ingredient that distinguishes the bosonic and 
fermionic MERA approaches. Remarkably, the cost of 
simulations does not depend on the particle statistics, 
but only on the amount of entanglement in the system. 

Sec. |IV| presents benchmark results. First a small lat- 
tice made of 6 x 6 sites is analyzed with the fermionic 
MERA and a fermionic tree tensor network (TTN), to 
confirm that the ansatz can accurately represent ground 
states. Both free and interacting systems are analyzed. 
Then much larger lattices, with up to 162 x 162 sites and 
periodic boundary conditions, are addressed in order to 
demonstrate the scalability of the present approach. Fi- 
nally, Sec. |V| contains some conclusions, and the appen- 
dices A, B and C provide some additional details. 

The present approach to account for fermions in a ten- 
sor network is equivalent to the one presented in Rcf. 
|2"51 We comment on this equivalence in appendix [D] The 
present form of the approach, however, makes its gen- 
eralization to other tensor network algorithms, such as 
PEPS (see also Ref. l29)l straightforward, as illustrated 
in Ref. 30 



II. BOSONIC MERA REVISITED 



A. Quantum circuit and renormalization group 
transformation 



whose output wires correspond to the sites of the lat- 
tice £o as depicted in Fig. [I] We first focus on the 
ternary ID MERA scheme introduced in Ref. [TSJ and 
then on its generalization to the 2D case. Several uni- 
tary gates transform the untentangled state 10)®^ into 
a state \^f) e V£ Q . We distinguish between two types 
of gates, isometries w and disentanglers u, each only in- 
volving a small number of input and output wires. A 
disentangler u is a map 
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with iy®2 the identity operator in V® , and an isometry 
w is a map 



with Jy the identity operator in V. 
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FIG. 2: (Color online) Disentanglers (a) and isometries (b) 
in the MERA are isometric gates, c) We usually draw the 
isometry without the incoming wires that are fixed to |0). 



Consider a lattice £q of N sites, where each site is de- 
scribed by a local Hilbcrt space Vo of finite dimension d. 
The MERA is an ansatz to describe certain states \^f) of 
the total Hilbert space Yc — Yf N , such as the ground 
state of a local Hamiltonian. The ansatz is efficient in 
the sense that the number of parameters required to en- 
code a state of a translation invariant system is only of 
order 0(log(N)x q ), with \ a refinement parameter (see 
below) and q a small integer number. This is in contrast 
to the dimension d N of the Hilbert space which grows 
exponentially with N. 




FIG. 1: (Color online) The ID MERA represented as a quan- 
tum circuit. It consists of different types of isometric gates: 
disentanglers (squares) and isometries (triangles). From the 
perspective of a renormalization group transformation the lat- 
tice £t_i is mapped to a coarse-grained lattice C T by applying 
a layer of disentanglers and isometries (cf. Fig. [3|. 



The MERA can be regarded as a quantum circuit 



From the perspective of a renormalization group trans- 
formation (see Ref. [17}, an isometry coarse-grains three 
sites into one effective site. A disentangler u acts across 
the boundary of two blocks of sites to reduce the amount 
of short-range entanglement between the blocks.^ The 
application of a layer of disentanglers followed by a layer 
of isometries describes a mapping of the lattice £ T _i into 
a coarse-grained lattice C r , as shown in Fig. [3| The local 
dimension of each coarse-grained site is denoted by x- 
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FIG. 3: (Color online) The real-space renormalization group 
transformation of the ID MERA (ternary scheme). 



Another key feature of the MERA is its causal struc- 
ture. The past causal cone of an outgoing wire s at time t 
is defined as the set of gates and wires that can affect the 
state in (s,t). A MERA is a quantum circuit for which 
the past causal cone of any location (s,i) in the circuit 
has a bounded width, i.e. involves only a constant number 
(independent of N) of wires at any previous time t' < t. 
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FIG. 4: (Color online) Top: The causal cone of an operator 
(oval) given by the shaded area involves only a small number 
of gates. Bottom: All gates outside the causal cone annihilate 
and we are left with a much simpler circuit. 



These key properties of the MERA enable an efficient 
calculation of expectation values of local observablcs, 
( , F|0| , I'), since only gates included in the causal cone 
of the operator O (i.e. the causal cones of the wires con- 
nected to O) have to be taken into account. All other 
gates can be replaced by identity operators thanks to 
Eqs. Q and as illustrated in Fig. [I] 

B. Superoperators and environments 

It has been shown in Ref. [THlthat the systematic eval- 
uation and manipulation of a MERA boils down to the 
calculation of several small diagrams which fall into three 
classes: ascending superoperators, descending superoper- 
ators, and environments. All these diagrams can be con- 
structed from the three generating diagrams shown in 
Fig. [5^)-c), as we explain in the following. 

Ascending superoperators - An ascending superoper- 
ator A transforms a two-site operator O t _i defined on 
lattice £ T -i into a two-site operator T on the coarse- 
grained lattice C T , as shown in Fig. [5ji). There are three 
structurally different ascending superoperators, which re- 
sult from the three generating diagrams in Fig. pa)-c) 




FIG. 5: (Color online) a)-c) The three generating diagrams 
of the ID MERA. d) An ascending diagram resulting from c) 
by erasing the two-site density matrix p T . e) A descending 
diagram obtained from c) by taking away the operator O t ~i- 
f) An environment for the disentangler u created by erasing 
the (upper) disentangler from c). An environment for w can 
be obtained in a similar way. 



by taking away the density matrix p T from each dia- 
gram. Repeated application of the corresponding ascend- 
ing superoperator to 0$ creates a sequence of increasingly 
coarse-grained operators {Oo,Oi, ...,Ot}- 

Descending superoperators - A descending superop- 
erator T> maps a two-site density matrix p T on the lat- 
tice C T into a two-site density matrix on the finer (less 
coarse-grained) lattice L T -i, as illustrated in Fig. [5^). 
The three different descending superoperators can be ob- 
tained from the generating diagrams by erasing the op- 
erator T _i from each diagram. Repeated application 
of the corresponding descending superoperator to pr-i 
creates a sequence of increasingly finer two-site density 
matrices {px-i, Pt-2, — j Po}- Note that px—i is obtained 
by joining the top-isometry with its conjugate. 

Environments - There are several ways to optimize 
the MERA. In this work we used the algorithm from Ref. 
1181 which is based on iterative optimization of individual 
gates. The optimization of, e.g., a disentangler involves 
calculating its three different environments, which are ob- 
tained by erasing the (upper) disentangler from the gen- 
erating diagrams, as shown for example in Fig. [5J 7 ). Note 
that an isometry has six different environments, three re- 
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suiting from erasing the left isometry from the generating 
diagrams, and the other three from erasing the one on the 
right (see Ref. [151 for more details). 

C. Diagrams 

Once we determined the diagrammatic representation 
of the ascending/descending superoperators and environ- 
ments we can evaluate the diagram by contracting the 
corresponding tensor network. We start by identifying 
all elements appearing in a diagram, which are summa- 
rized in Fig. [6j 

1. Elements in a diagram 



degrees of freedom. In general, we replace each crossing 
by a swap gate B which accounts for the exchange pro- 
cess (see appendix [A]) . In the bosonic MERA this gate 
is simply the identity, i.e. B % £? = 5i 1 j 1 5i 2 j 2 , because 
the bosonic wavefunction is symmetric under exchange 
of two particles. As a consequence the crossings can sim- 
ply be ignored. However, this will no longer hold in the 
fermionic MERA as we will see Sec. IIIII 

Lines - A single line in a diagram corresponds to the 
identity 7y of the vector space V. A line connecting two 
tensors describes how they are multiplied together, when 
the diagram is contracted, as we explain next. 

2. Contraction of a diagram 
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FIG. 6: (Color online) Elements in a diagram of the ID 
MERA. Each shape represents a tensor: a) isometry, b) dis- 
entangler, c) Hamiltonian, d) density matrix, e) A crossing 
of lines corresponds to a two-body gate (which is simply the 
identity in the bosonic MERA). f) A single line corresponds 
to the identity Iy of the vector space V. 



Shapes - Each shape represents a tensor (a multidi- 
mensional array) with a rank equal to the number of legs. 
The entries of a tensor are given by the expansion coeffi- 
cients of the corresponding gate (or Hamiltonian/dcnsity 
matrix) in the local bases of its legs. For example, a 
general two-body Hamiltonian term can be expanded as 



H 



E 



(3) 



where each sum goes over all basis states of the local 
Hilbert space of each leg. The four-leg tensor associated 
to H is H* 1 ] 2 . 

3l]2 

Line crossings - A diagram may involve line crossings, 
e.g. the ID MERA in the case of periodic boundary con- 
ditions (or for a Hamiltonian with next-nearest neighbor 



interaction), or the 2D MERA as we will sec in Sec. II D 



Each crossing corresponds to an exchange of the degrees 
of freedom carried by the individual lines. The implica- 
tions of this exchange depend on the statistics of the basic 



FIG. 7: (Color online) a) Multiplication of two tensors H and 
u on the legs connected by the line labelled by i%. b) A line 
connecting two legs of the same tensor corresponds to a trace 
(see text). 



A tensor network is contracted by a sequence of pair- 
wise multiplication of tensors. Two tensors are multiplied 
together according to the lines connecting the legs of the 
tensors. For example, the multiplication in Fig. [7^) of 

h 2 h 3 
"3132 J _i 2 i3 

tensor Hu given by 



-Hi 1 !? with u" 2 ," 3 on the leg labelled by i^ leads to a new 



\Hu\ ilh2fl3 = V* H h 

f" U l 3132*3 31 



t2„h 2 h 3 

32 1213 



(4) 



A special case is the trace, which is represented as a line 
connecting two legs of the same tensor, as for example 



3233 ' 



(5) 



illustrated in Fig. jjp). An example of a full contraction 
of a tensor network is shown in Fig. [8] 

Note that the computational cost depends on the order 
in which the pairwise multiplications are implemented. 
In a practical implementation it is therefore crucial to 
determine the order which minimizes the computational 
cost (and/or memory requirements). The computational 
cost to multiply two tensors A and B connected by l c 
legs, is given by x? a+1b ~ 1c , where Ia (Ib) is the number 
of legs of tensor A (B), and we assumed that each leg has 
the same dimension \- The scaling of a MERA algorithm 
is dominated by the largest cost in the contraction of a 
diagram. The cost in memory scales with , with 
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FIG. 8: (Color online) Contraction of the tensor network in 
Fig. |HJi) by a sequence of pairwise multiplication of tensors. 



Imax the tensor with the biggest number of legs occur- 
ring during the contraction. For the ID ternary MERA 
the computational cost scales as 0(\ 8 ), and the cost in 
memory as 0(x 6 )P^ 
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FIG. 9: (Color online) The real-space renormalization group 
transformation of the 2D MERA. 



D. 2D MERA 

There are s everal ways to realize a MERA in 
2D j i4 | iM | 2i | 23 |HlHere we focus on a "9-to-l scheme" where 
a site of C T corresponds to a block of 3 x 3 = 9 sites of 
£ T -i, and with disentanglers that do not overlap, as de- 
picted in Fig. M The computational cost of this scheme 
scales as 0(xji and the cost in memory as 0(x 12 )- 

Conceptually, one proceeds in the same way as in the 
ID MERA, i.e. one determines ascending/descending su- 
peroperators and environments. The ascending superop- 
erator maps a 4-body plaquette operator O r _i on the 
lattice £ T _i into a plaquette operator O t on the lattice 
C T . But a 2-body operator, depending on its location 
in the lattice, may be mapped into a 4-body operator 
on a higher level. We therefore focus here on plaquette 
operators. (Note that in case of a two-body Hamilto- 
nian we can either treat the lowest layer differently than 
the higher layers, or express the hamiltonian as a sum of 
plaquette operators from the start). 

Figure [10] shows one particular generating diagram, 
from which we can obtain ascending/descending super- 
operators or an environment by eliminating the corre- 
sponding tensor, as explained for the ID case. There are 
9 different generating diagrams, corresponding to the 9 
different positions of the Hamiltonian with respect to the 
basic 3x3 block. 

Note that the basic diagrams of a 2D MERA arc 2 + 1 



dimensional objects. In Figure 10 we chose one particular 
way to map the diagram onto a plane, i.e. 1+1 dimen- 
sions like the ID MERA. This mapping is not unique, 
but for any choice, some of the lines in the diagram cross 



each other. As already mentioned, these crossings can 
be ignored in the bosonic case. However, they will play 
an important role for the fermionic MERA, as we will 
explain in the next section. 



III. FERMIONIC MERA 

The essential difference between a fermionic and a 
bosonic system lies in the symmetry of the wavefunction 
under the exchange of two particles. Exchanging two 
bosons leaves the wavefunction invariant, whereas when 
exchanging two fermions the wavefunction is multiplied 
by — 1. More generally, exchanging an odd number of 
fermions living on a (coarse-grained) site i' with an odd 
number of fermions on j' leads to a negative sign. 

All basic concepts introduced for the bosonic MERA 
still hold for the fermionic MERA, i.e. the gates are iso- 
metric and the causal cone is the same as in the bosonic 
MERA (see appendix [C). All we need to do is to use par- 
ity preserving tensors, and introduce a fermionic swap 
gate, which implements the fermionic exchange proper- 
ties, as we explain in the following. 



A. Z2 symmetry 

A property of any fermionic Hamiltonian H (and more 
generally any fermionic observable) is that it preserves 
parity, i.e. [P,H] — 0, with P — (— 1) N the total parity 
operator, where N measures the total number of parti- 
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a and b into one site c can be split into a fusion of the two 
sites (blocking) followed by a truncation of the combi ned 



Hilbert space V — V a <8> V^, as illustrated in Fig. 11 
The fusion rules describe how the individual sectors of 
V = V^ + ' © result from combining the sectors of V a 
and Vfc: 

= (YWgyW)®^-)®^-'), (7) 

v<-) = (v^ ) ®v{ i +) )e(vi +) ®v^ ) ). (8) 



Note that the truncation is performed separately in each 

Finally, a fusion 



sector, V (+ ) -> ¥ { c +} and 
of N sites can be decomposed into N 
as illustrated in Fig. [llj>). 
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FIG. 11: (Color online) a) An isometry that coarse-grains two 
sites a and b into one site c can be decomposed into a fusion of 
the two sites followed by a truncation of the combined Hilbert 
space V. b) An isometry that coarse-grains three sites into 
one can be decomposed into two subsequent fusions, followed 
by a truncation. 



FIG. 10: (Color online) Example of a generating diagram 
of the 2D MERA projected onto 1+1 dimensions. On the 
bottom we define the correspondence between the gates in 
this figure with the ones from Fig. [9] as indicated by the 
numbers. The picture on the bottom left shows the location 
of the 4-body plaquette Hamiltonian (oval) in the lattice. 



We emphasize that also many bosonic systems ex- 
hibit a Z2 symmetry, which can be incorporated into the 
bosonic MERA in the same wayP*l In general, exploiting 
symmetries increases the efficiency of a simulation. How- 
ever, for fermions, the parity symmetry also plays an im- 
portant role for the implementation of the fcrmionic swap 
gate which we introduce in the next section. 



cles in the system. This Z2 symmetry stems from the 
fact that fermions can only be created or annihilated in 
pairs. We incorporate this symmetry into the MERA by 
enforcing all tensors to be parity preserving. A tensor 
r Ti 1 i- 2 ...i M preserves parity if 

T ili2 ... iM = 0, \iP{i x )P{i 2 )...P{i M )^l, (6) 

where P(ife) € {—1, 1} denotes the parity of the state la- 
belled by ik. The local Hilbert space of a (coarse-grained) 
site is decomposed into a space with even parity (+) and 
one with odd parity (-), i.e. V = V (+) ® V (_) . Each 
basis state in V is labeled now by a composite index 
j = (p,a p ), where p 6 {+, — } specifies the parity sec- 
tor and aP enumerates the states in the subspace V^. 
This decomposition allows us to identify the parity of a 
state very easily, and it also leads to a block structure of 
the tensors (similarly to a block diagonal matrix). 

Fusion rules - An isometry that coarse-grains two sites 



B. Fermionic swap gate 
parity (+) (+) (+) (-) (-) (+) (-) (-) 




parity (+) (+) (-) (+) (+) (-) (-) (-) 
sign +1 +1 +1 -1 

FIG. 12: (Color online) The fermionic swap gate implements 
an exchange of fermions. Exchanging an odd number of 
fermions on one site with an odd number of fermions on an- 
other site leads to a negative sign factor. 
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As explained in Sec. II C two crossing lines i and j 



in a tensor network are nothing but a graphical repre- 
sentation of an exchange process (or a swapping). As a 
consequence of the antisymmetry of the fcrmionic wave- 
function, a prefactor of —1 appears if both lines carry a 
state with odd parity (odd number of particles), as illus- 
trated in Fig. [12] We replace each crossing by a gate 
B that accounts for this exchange process (see appendix 
lAl): 



with 



(9) 



(10) 



only depending on the parities of the states i\ and i 2 - 
The function S evaluates to —1 if both parities are odd, 
and +1 otherwise. 

Having parity preserving tensors allows us to take a 
line and "jump" over another tensor, as illustrated in Fig. 
[T3| We demonstrate the validity of this transformation 
in appendix[B] Before contracting the tensor network, we 
rearrange the lines and tensors in such a way that each of 
the resulting fermionic swap g ates ca n be absorbed into a 

M The resulting tensor 
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single tensor, as shown in Fig. 
network can then be contracted in the same way as in 
the bosonic case. Note that the computational cost of 
absorbing a fermionic swap gate into a tensor with I legs 
is only of order x ■ Therefore, this cost is subleading and 
the overall cost of the algorithm is essentially the same 
as in the bosonic case. 

In other words, we map a non-planar tensor network 
to a planar one (i.e. a network without line crossings) by 
replacing the crossings by fermionic swap gates. We can 
modify the resulting planar network by "jump" moves in 
such a way that the resulting fermionic swap gates do 
not increase the leading cost of a contraction, compared 
to the bosonic case. 

The fermionic MERA presented in this paper may look 
different than the one introduced in Ref. HH1 which is 
based on the Jordan- Wigner transformation to map the 
fermionic system into a bosonic one. However, it is im- 
portant to point out that the two approaches describe 
the same MERA (see appendix |A| . 



IV. RESULTS 

In this section we present benchmark results for the 
fermionic 2D MERA, and also for the 2D tree tensor 
network (TTN) which corresponds to the 2D MERA 
without disentanglers. Thanks to its simpler structure, 
a larger value of \ is affordable. However, in contrast to 
the MERA the TTN is not scalable, i.e. in general \ nas 
to increase exponentially with system size to account for 
the accumulation of short-range entanglement across the 
boundary of a block. 






I absorb 



FIG. 13: (Color online) An example diagram involving line 
crossings. Thanks to the parity symmetry of each tensor we 
are allowed to "jump" with a line over a tensor in order to 
simplify the tensor network. Lines are moved around in such 
a way that each fermionic swap gate can be absorbed into a 
single tensor. The resulting tensor network is contracted as 
in the bosonic case. 



We first consider an exactly solvable model of non- 
interacting spinless fermions in two dimensions given by 
the Hamiltonian 



(rs) 



+4c r -7(44+c s c 



]-2A J2 4<V 



(11) 



with A the chemical potential and 7 the pairing potential. 
The phase diagram of this model (see Fig. 14) exhibits 
a critical (p-wave) superconducting phase for 7 > 0, 
< A < 2 with two gapless modes, and a gapped su- 
perconducting phase for 7 > 0, A > 2) 38 * 39 l For 7 = the 
model corresponds to a free fermion system, i.e. a metal 
(with a one dimensional Fermi surface) for < A < 2 
and a band-insulator for A > 2. A metallic phase is also 



found for 7 > and A = 0. The lower panels in Fig. 14 
present the error in the ground state energy of a 6 x 6 
system as a function of 7 and A, for increasing values of 
X- Both TTN and MERA reproduce several significant 
digits of the exact solution. The middle panels show the 
entanglement entropy of half the system, 



5*1/2 — — ^ 

k 



log 2 A fc , 



(12) 



with A/c the eigenvalues of the reduced density matrix of 
half the system. The accuracy of the energy is clearly cor- 
related with the amount of entanglement in the system, 
i.e. the accuracy decreases with increasing Si/ 2 - Ac- 
curate results are also obtained for correlators, C(r) = 



4 

To 1 



as shown in Fig. 15 Note that the 6x6 sys- 



tem corresponds to a MERA with only one single layer of 
isometries and disentanglers, which has a computational 
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FIG. 14: (Color online) Top left panel: Phase diagram of the 
free fermion model (111. Lower panels: Error in the ground 



state energy obtained from TTN and MERA simulations of 
a 6 x 6 lattice with periodic boundary conditions. The lines 
correspond to different values of the refinement parameter 
X, as indicated in the legend in the top right panel. The 
accuracy of the simulation results depends on the amount of 
entanglement in the system, which is measured here by the 
entanglement entropy of half the system, S1/2, plotted in the 
middle panels. 



cost that scales as 0(\ 4 ) (times a factor depending on 
the local dimension d in the lattice Co)- This allows us 
to use a larger value of \ than in the MERA with several 
layers (large systems, see below). 



Next we consider the same Hamiltonian (11) with an 
additional nearest-neighbor interaction, 



Hint — Hfre 



(rs) 



C r C r C a Cc 



(13) 



which can no longer be solved analytically. We empha- 
size that the algorithm does not require any particular 
modification in order to deal with the interaction, since 
an arbitrary 2-body Hamiltonian can be used as an input 
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FIG. 15: (Color online) Top panels: Correlation function 
C(r) = (c?, cpQ+p) for 7 = 1, A = 1.5 (left) and for 7 = 1, A = 
2.5 (right). The positions f are the same as indicated in the 
bottom plot. Lower panels: The difference between the sim- 
ulation result and the exact analytical solution for different 
values of the refinement parameter \. 



convergence of the energy with x f° r different interaction 
strengths V. For small (V « 1) and large (V >> 1) 
interaction we find a similar convergence behavior as in 
the non- interacting case. In both cases 5 , 1 / 2 i s relatively 
small, as shown in the upper panel of Fig. 16 For an in- 



teraction strength of the order of the hopping amplitude, 
V ~ t = 1, the convergence with x is slower but rs 4 dig- 
its of accuracy are still achieved for large x- Accordingly, 
the amount of entanglement in the system (measured by 
S1/2) is large in this parameter region. As another ex- 
ample of a correlation function we computed the pairing 
amplitude 



p(k) = <4ci t >, 6 



1 



E 



4exp(zkr). (14) 



to the simulation. The lower panels in Fig. 16 show the 



Figure [17] shows the total pairing amplitude Ptot — 
^ k |P(k)| as a function of V, for two sets of parame- 
ters for A and 7. Also this quantity converges to the 
exact solution with increasing x i n the exactly solvable 
case, V = 0, as shown in the inset. A similar convergence 
behavior is observed in the interacting case (not plotted) . 
A small interaction amplifies the total pairing amplitude, 
whereas a large interaction tends to suppress the pairing. 
The sudden jump of both curves around V rs 2.2 could 
indicate a first order phase transition. A weaker feature 
is found for V rs 1.5 (crosses) and V rs 1.2 (dots). 

Finally, we show that the MERA is scalable in two 
dimensions. Figure 18 shows the relative error of the 
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FIG. 16: (Color online) Convergence of the ground state en- 
ergy of interacting spinless fermions on a 6 x 6 lattice with 
periodic boundary conditions and y = 1, A = 2 for different 
interaction strengths V. The plot shows AE — E x —E min , the 
difference between the energy as a function of \ an d the best 
(lowest) energy obtained by the simulations. Open squares 
are obtained by the TTN, filled circles by the MERA. 
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FIG. 18: (Color online) The relative error of the energy as 
a function of system size obtained by the 2D MERA is of 
the same order of magnitude for small systems as for large 
systems, even in the critical regime A < 2. The simulations 
are done for x — 4 up to a system size 162 x 162, and V = 0. 
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FIG. 17: (Color online) Total pairing amplitude Ptot as a 
function of the interaction strength V of the spinless fermion 
model (13), obtained by a TTN with \ = 200 - 300. Larger 
values of x produce corrections between 10~ 3 (for V ~ 1.5) 
and 10 -4 . The inset shows the convergence of Ptot with x m 
the exactly solvable case, V = 0. 



energy as a function of system size up to 162 x 162 for 
the non-interacting case V = For a fixed x — 4 the 
relative error is of the same order of magnitude for small 
systems as for large systems, even in the critical regime 
A < 2. The system size can easily be increased by adding 
more layers of the MERA, with a cost that only grows 
logarithmically with the system size for a translational 



invariant system 
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To summarize, we have shown that fermionic systems 
can be addressed in a very similar way as bosonic systems 
within the formalism of entanglement renormalization. 
We explained how to modify the bosonic MERA in order 
to deal with fermionic degrees of freedom. To do this 
we incorporate a Z2 symmetry into the MERA by using 
parity preserving tensors, and introduce a fermionic swap 
gate to account for the exchange of fermionic degrees of 
freedom, whenever two lines in the tensor network cross. 
We showed that a fermonic tensor network can be trans- 
formed in such a way that the fermionic swap gates do 
not increase the complexity of a contraction. Thus, an 
important result is that the complexity of the fermionic 
MERA is the same as for the bosonic MERA. The present 
formalism to deal with fermionic systems was originally 
developed specifically for the TTN and the MERA in 
mind but, in its present formulation, can be applied to 
arbitrary tensor networks. The steps to be followed are 
surprisingly simple: given a tensor network ansatz, such 
as PEPS, one must choose Z2 symmetric (i.e. parity 
invariant) tensors and, when contracting the tensor net- 
work, replace crossings with fermionic swap gates. This 
procedure is exemplified in Ref. [30] for infinite PEPS. 

Here we have presented benchmark results for the 2D 
MERA and the TTN for spinless fermions, both non- 
interacting (exactly solvable) and interacting. We have 
also shown that the 2D MERA is scalable by simulating 
lattices made of up to 162x162 sites. 

In general, the efficiency of the MERA depends on 
the amount of entanglement in the ground state of the 
system. Accordingly, as discussed in Ref. [U for free 
fermions, gapped systems appear typically as the eas- 
iest to simulate. They are followed by critical phases 
with a finite number of zero modes (e.g. Dirac modes), 
which are more entangled but still follow an area law 
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for the entanglement entropy,^ 2 | 33 | 34 | 3S | which a MERA 
with the same x a t each level of coarse-graining can 
reproduce.^ The most challenging systems are metals 
with a one dimensional Fermi surface, i.e. an infinite 
number of zero modes. These are the most entangled 
systems, wit h a mu ltiplicative logarithmic correction to 
the area law F ' 33 ' 34 ! 3 * ! 

The efficiency of the algorithm can be substantially 
improved by making use of symmetries (e.g. SU(2), 
U(l), etc.) of a mod elPj] a nd by variational Monte Carlo 
sampling techniques! 41 * 42 ! This is important in order to 
increase the maximal affordable x, which for large sys- 
tems is still small at present. A higher accuracy can also 
be achieved by choosing an optimal structure of the 2D 
MERA depending on the problem considered. An exam- 
ple of an improved coarse-graining scheme was presented 
in Ref. [24] 

We believe that the fermionic MERA will help to 
shed new light into long standing questions in strongly 
correlated fermion systems. Work in progress includes 
the study of the ground state phase diagram of the tJ 
and the Hubbard model, and generalizations to anyonic 
systems. S3 

We thank R. Pfeifer and L. Tagliacozzo for useful dis- 
cussions, and S. Haas, L. Ding and N. Ali for clarifi- 
cations concerning the free fermion model (111. Support 
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NOTE: Short after Ref. |28| (of which the present pa- 
per is an extended version) was made available online, a 
largely equivalent approach has been independently pre- 
sented by C. Pineda, T. Barthel and J. Eisert in Ref. 



APPENDIX A: THE FERMIONIC SWAP GATE 

In this appendix we show that the fermionic swap gate 
implements the anticommutation of fermionic operators, 
and make connection to the Jordan- Wigner transforma- 
tion used in Ref. [2U 

In a fermionic MERA the wires in the quantum cir- 
cuit carry fermionic degrees of freedom, and all gates 
and operators can be expanded as a product of fermionic 
creation and annihilation operators. These operators an- 
ticommute instead of commuting as in the bosonic case. 
To illustrate this essential difference between a fermionic 
and a bosonic MERA (for hardcore bosons) we consider 
the computation of a matrix element of an operator O 
acting on sites 1 and 3 in a three-site system, shown in 
Fig. [19] The full Hilbert space is spanned by the basis 
states 



Stil 't*2 -t*3 I 



where ik € {0, 1} and fit creates a particle on site k, 
obeying the commutation relations 



[Ci,C 



3l± 



0, 



At 




ii44lo) 



±ix44lo> 

±tc\c 3 I 1 cl&l\0) = 
= ±tc\l 3 cl\0) 

X 

±fcfeJ,l3|0) 
110|(±i)|110) =±t 



FIG. 19: (Color online) A diagram representing the matrix 
element {ip' \0\tf)) . Reading the diagram from top to bottom, 
we obtain a prescription of how to calculate the matrix ele- 
ment. Crossing lines imply that the operators carried by the 
line are exchanged, which results in a negative sign in the 
fermionic case if two creation (or annihilation) operators are 
exchanged. The resulting matrix element is +t in the bosonic 
and —t in the fermionic case. 



In the bosonic case (+) operators on different sites com- 
mute, whereas in the fermionic case (— ) they anticom- 
mute. Let us consider the example of a hopping operator 
6 = tc\c 3 , and the states = |011) and (ij>'\ = (110|. 
Figure [19] provides a graphical prescription of how to 
compute the matrix clement, 



{i>'\0\i>) = (110|tc{c 3 |011) = i(0|c 2 c 1 c t 1 c 3 cjcj|0). (A3) 

For a bosonic system this matrix element is simply t, be- 
cause operators on different sites commute. However, in 
the fermionic case a minus si gn a ppears from exchanging 



with c\, as shown in Fig. 19 Thus, a line crossing is 



a graphical representation of the exchange of operators, 
and a negative sign results if both lines carry an odd 
number of fermionic creation/annihilation operators. 

More generally, consider an arbitrary (parity preserv- 
ing operator) O acting on sites 1 and 3 and a state 
expanded in the local basis, 



The expectation value is given by 



\0H) = 



(A4) 



£ 

H*2*3 
31J3 

x(0|(c 2 ; 3 c-ci 1 )(c^c^C3 3 ci 1 )(cr i 4-c 3 -)|0) 



(Al) In the bosonic case the second line in Eq. (A4| is al 



{1,2,3}. (A2) 



ways 1, because the bosonic operators commute, and the 
expected value is simply obtained by multiplying the ten- 
sors in the first line together. In the fermionic case we 
again have to swap operators as indicated by the cross- 
ing lines in Fig. 20 1) . More conveniently, we replace each 
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a) 



b) 



^Jl«2j3 



FIG. 20: (Color online) a) A general expectation value of a 
two-body operator acting on sites 1 and 3. As usual, the 
crossing of lines implies exchanging the operators carried by 
the lines, b) Mapping to a tensor network by replacing each 
line crossing by a swap gate, which implements the anticom- 
mutation of fermionic operators. 





MERA presented in this paper. A third equivalent ap- 
proach (but yet another point of view) is to change the 
Jordan- Wigner order of the lattice sites such that op- 
erator O acts on contiguous sites, as illustrated in Fig. 



21] The gate S to change the Jordan- Wigner order cor 
responds to the swap gate B introduced in Eq. (10 1 



except that the lines do not cross, i.e. . ... , , 



We 



summarize the three equivalent approaches in appendix 



APPENDIX B: PROOF OF THE "JUMP" MOVE 



The "jump" move introduced in Sec. IIIB allows us 
to drag a line over a tensor, as for example shown in Fig. 
p~3| The validity of this rule originates from the particular 
form of the swap gate and from the fact that all tensors 
preserve parity, as we explain in the following. 



crossing by a swap gate introduced in Eq. (10 1, which 
accounts for the anticommutation rules of the fermionic 
operators, so that the expectation value becomes 

(1>\d\4>) = £ r il233 Bl:;iO^B^ Hl2l3 . (A5) 



31J3 



Therefore, replacing each line crossing by a swap gate 
transforms the diagram in Fig. [20^ ) into the tensor net- 
work shown in Fig. 20 :>), which we can contract as ex- 
plained in Sec. |II C| 



at 



b) 



C) 



O 





<B 








FIG. 21: (Color online) Three equivalent approaches to com- 
pute an expectation value of a fermionic two-body operator 
acting on sites 1 and 3: a) by using swap gates as presented 
in this work, b) by using a Jordan- Wigner transformation of 
the fermionic operator O into a three-site spin (bosonic) op- 
erator O, c) by changing the Jordan- Wigner order of sites 2 
and 3 by a gate 5* (cf. text), so that the operator O acts on 
contiguous sites (with respect to the Jordan- Wigner order as 
indicated by the numbers). 



If we incorporate the swap gates into the operator O 
we end up with the three-site operator shown in Fig. 
20 }) , which is nothing but the Jordan- Wigner transfor- 
mation of operator O, i.e. the fermionic operator mapped 
to (bosonic) spin variables. This approach was used to 
introduce the fermionic MERA in Ref. [551 but it is im- 
portant to point out, that it is equivalent to the fermionic 




d) i 




e) h 



f) h k 




FIG. 22: (Color online) a) A parity preserving tensor remains 
invariant when applying the parity operator P on each leg. 
b) This equality follows from P 2 = 1, i.e. the parity operator 
squared is the identity, c) Two subsequent line crossings is 
equivalent to the identity, d) and e) "Jump" move over a 
tensor with two legs, corresponding to Eq. (Bll. f) and g) 



"Jump" move over a tensor with three legs, explained in Eq. 
(B2I. 



First note that, by definition, a parity preserving ten- 
sor remains invariant when the parity operator P acts on 
all of its legs, as illustrated in Fig. 22 i). It then follows 
from P 2 = 1 that acting with P on a set of legs of a ten- 
sor is equivalent to acting with P on the complementary 
set of legs of the tensor, see Fig. 22 }). Next we consider 



a two-legged tensor A and a swap gate B (see Eq. (lOl) 
shown in Fig. |22]i). Contracting A and B amounts to 



*2 

= S jujl A^(P(h),P(ii)) = E B S;4 (Bl) 



where we made use of the parity symmetry of A, i.e. that 
acting with the parity operator P on the lower leg of A 
is equivalent to acting with P on the upper leg. The 
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final expression in Eq. (Bl I is nothing but the swap gate 



now acting on the upper leg of tensor A as shown in Fig. 
22:). In a similar way one can proof the "jump" move 
for a three-legged tensor shown in Fig. 22 ') and g) by 
contracting the tensor networks and using 

S(P(h),P(i 3 )) = 5(P(j 1 ),P(i 1 )P(z 2 )) (B2) 
= S(P(j 1 ),P(i 1 ))5(P(i 1 ),P(i 2 )). 

The first equality again follows from the parity symme- 
try, and the second equality can be easily verified. This 
property can be extended to tensors with more than 3 
legs by applying the last identity recursively, i.e. 

S(P( J1 ),Y[P( lk )) = Y[S(P( J1 ),P( lk )) (B3) 

k k 

Finally, we show that two subsequent line crossings of the 
same lines is simply the identity, as shown in Fig. 22;). 



This follows from multiplying two swap gates together, 

(B4) 



1 k 1 k 2 - hh feifos 



= S ilM S l2M (S(P(n),P(t 2 )f 

$i 1 ,ki$i2,k2 ■ 

By combining above rules an arbitrary jump move can 
be performed. 



APPENDIX C: CAUSAL CONE OF THE 
FERMIONIC MERA 

It is important to notice that the fermionic MERA has 
the same causal cone as the bosonic MERA (cf. Fig. [4]). 
This can be understood as a consequence of the "jump" 



move, as shown in Fig. 13 An arbitrary diagram can be 
modified in such a way that no line crosses the outgoing 
wires of a gate that lies outside the causal cone. The gate 
with its conjugate can then be replaced by the identity 
thanks to Eqs. dTl and ([2]), as done for example in the 
diagram in Fig. [13] with the isometry in the middle. 

Another way to arrive to this conclusion is to consider 
the MERA where we expand each gate and the Hamilto- 
nian in its local fermionic basis. An expectation value of 
an operator O is of the form 



(d)=(o\wl u 



rjv„' 



.u\Ou\. . .un u wi. . -Wn^ |0) (CI) 



where we consider only one layer of the MERA for 
simplicity, and N w and N u are the number of isometries 
and disentanglers in the layer, respectively. Because 
each gate is parity preserving, its expansion consists only 
of terms with an even number of creation/annihilation 
operators. As a consequence two parity preserving gates 
with disjoint supports (i.e. gates acting on different 
sites) commute. Thus, to simplify Eq. |C1| we can 
first commute each disentangler u^. that lies outside 
the causal cone of operator O to the left to annihilate 



with its conjugate, and then proceed similarly with the 
isometries Wk, so that only gates inside the causal cone 
of operator O remain. For example, if Hi lies outside 
the causal cone, [O, iii) = and ii\u\ = 1 annihilate. 
This generalizes straightforwardly to several layers of 
isometries and disentanglers as in Fig. [4] 



APPENDIX D: EQUIVALENT APPROACHES 

In this appendix we briefly review the formulation of 
the fermionic MERA originally presented in Ref. dHl to- 
gether with an alternative formulation also outlined in 
Ref. 1281 and compare it to the simplified formulation 
presented in this paper. 

1 ) Fixed Jordan- Wigner order: The formalism intro- 
duced in Ref. HHHs based on the Jordan- Wigner transfor- 
mation (with a fixed Jordan- Wigner order), which maps 
all fermionic operators into spin operators with string of 
Z's (the a z Pauli matrix in case of a two dimensional local 
Hilbert space). An example of such an operator is shown 
in Fig. 



21 d). These strings of Z's are coarse-grained 



locally by using fermionic disentanglers and isometries. 
The fermionic trace allows us to dispense with the string 
of Z's, as noticed in Ref. [28] For example, the reduced 
density matrix of sites 1 and 3 of a system with three 
sites is computed as 

P13 = ftr 2 (pi23) 



a=0,l 



n a \x^ Pl23 z 2 )v[ 



(1-eO 



(Dl) 



a=0,l 



where Vr°^ and Vr project onto the even and odd parity 
sectors of site r, and p\ 23 is the density matrix of the full 
system. Indeed, one can check that 



tr(A 1 Z 2 B 3 p 123 ) 
ti(C 1 I 2 D 3 p 1 2 3 ) 



tr(AiP 3( oi 3 ) 
tr(dD 3Pl3 ) 



(D2) 
(D3) 



for A, B parity changing operators (odd parity opera- 
tors) and C, D parity preserving operators (even parity 
operators) 

2) Changing the Jordan-Wigner order: The use of 
a fermionic trace amounts to effectively changing the 
Jordan-Wigner order in such a way that, in the new or- 
der, the sites to be kept after tracing out are contiguous 
sites. Similarly, before applying a specific operator, the 
Jordan-Wigner order is changed accordingly so that the 
operator acts on contiguous sites in the new order, as 
shown, for instance, in Fig. 21:). The same holds for 
isometries and disentanglers. This was pointed out in 
Ref. US] (see also Refs. 1341351) . 

3) Crossing of lines carrying fermionic degrees of free- 
dom: As discussed in appendix [X] the above approaches 
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are also equivalent to the one presented in this paper. 
The latter has the advantage of completely dispensing 
with the Jordan- Wigncr order, and its implementation 
on generic tensor network algorithms for ID lattices with 



periodic boundary conditions (e.g. MPS, TTN, MERA), 
and 2D lattices (e.g. PEPS, TTN, MERA) is straight- 
forward. 
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